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Final  Report 

DARPA  Grant  F49620-87-C-0065 
December,  1990 

General  Survey 

The  principle  goal  of  the  work  done  on  DARPA  Grant  No.  F49620-87-C-0065  has  been  to 
produce  new  algorithms  for  solving  problems  in  the  properties  of  fluids  and  solids  and  to  find 
new  techniques  for  the  parallelization  of  algorithms.  Foci  have  been  on  algorithms  for  the  non¬ 
destructive  testing  of  objects,  the  optimization  of  structures,  the  computational  analysis  of  phase 
change,  the  parallelization  of  linear  algebra  and  other  basic  algorithms.  Other  topics  have  been 
involved  because  of  their  intimate  connection  either  to  parellelizaton  (e.g.  queueing  theory),  to 
materials,  or  to  numerical  methods  for  relevant  p.d.e.’s  including  nonlinear  hyperbolic  methods 
themselves. 

This  report  consists  of  a  summary  of  the  topics  studied,  the  investigators  and  short  outlines 
of  their  achievements.  These  are  grouped  under  the  titles.  Algorithms  for  Materials,  Paralleliza¬ 
tion  of  Algorithms,  Other  topics.  A  general  description  is  followed  by  longer  reports  covering 
the  past  year  and  a  half  of  activity,  a  bibliography  of  work  published  and  a  list  of  researchers  and 
students  supported. 

Algorithms  for  Materials 

This  work  was  conducted  mainly  under  the  guidance  of  R.  Kohn,  J.  Goodman  and  C.S. 
Morawetz.  It  involved  as  well  G.  Allaire,  S.  Cox,  A.  Szepassy,  J.  Sylvester,  John  Strain,  M. 
Landman  and  Paula  Whitlock.  The  areas  were  structural  optimization,  electrical  impedance 
tomography. 

Structural  Optimization 

The  goal  of  structural  optimization  is  to  choose  the  geometry  or  composition  of  a  load- 
bearing  structure  so  as  to  optimize  its  performance  for  a  given  amount  of  material. 

S.  Cox  worked  with  M.  Overton  on  non-smooth  optimization  problems  arising  in  the 
optimal  design  of  columns  against  buckling.  Prior  work  had  been  based  primarily  on  necessary 
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conditions,  and  contained  mistakes  that  had  led  to  much  controversy  in  the  literature.  Cox  and 
Overton  corrected  those  mistakes,  and  showed  that  the  problem  can  actually  be  solved  by  a  direct 
optimization  scheme.  The  associated  mathematical  programming  problem  is  non-smooth, 
because  the  buckling  load  of  an  optimal  column  is  usually  a  multiple  eigenvalue  of  the  associated 
differential  equation.  This  work  was  reported  in  "On  the  optimal  design  of  columns  against 
buckling,"  by  S.  Cox  and  M.  Overton,  to  appear  in  SIAM  Journal  of  Mathematical  Analysis. 

R.  Kohn  and  G.  Allaire  explored  the  shape  optimization  of  two-dimensional  structures  in 
plane  stress.  They  implemented  a  new  approach,  based  on  the  introduction  of  "composite  materi¬ 
als"  as  design  elements.  This  amounts  to  permitting  subgrid  structure,  hence  to  enlarging  the 
apparent  design  space.  The  result  is  an  approach  to  structural  optimization  that  is  more  global  in 
character  -  with  fewer  local  minima  than  conventional  methods.  This  work  is  now  being 
prepared  for  publication  in  the  form  of  an  article  "Optimal  design  for  plane  stress  using  compo¬ 
site  materialst'•  by  G.  Allaire  and  R.  Kohn. 

R.  Kohn  also  explored  a  connection  between  optimal  design  and  coherent  phase  transitions. 
The  latter  is  actually  also  a  problem  of  optimal  design  -  in  which  the  quantity  to  be  optimized  is  a 
thermodynamic  free  energy.  The  "composite  materials"  discussed  above  correspond  to  the 
microstructures  of  phase  mixtures  created  by  coherent  phase  transitions.  Kohn  explored  the 
metallurgical  literature  in  this  area,  which  goes  back  to  the  1960’s,  and  explained  its  relation  to 
the  recent  work  on  structural  optimizaton.  This  work  is  reported  in  the  "The  relaxation  of  a 
double-well  energy,"  by  R.  Kohn,  to  appear  in  Continuum  Mechanics  &  Thermodynamics. 

Electrical  Impedance  Tomography 

The  goal  of  electrical  impedance  tomography  is  to  image  the  interior  conductivity  of  an 
object  by  menas  of  low-frequency  voltage  and  current-flux  measurements  at  the  boundary.  This 
is  a  promising  new  approach  to  nondestructive  testing.  It  has  also  been  proposed  for  a  variety  of 
biophysical  and  geophysical  applications. 

R.  Kohn  and  A.  McKenney  explored  the  implementation  of  a  new  variational  approach  to 
this  inverse  problem.  They  developed  a  Newton-type  minimization  scheme,  and  explored  the 
performance  of  the  method  in  a  two-dimensional  setting.  It  was  found  to  be  surprisingly  stable  in 
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the  presence  of  noise,  provided  that  the  boundary  currents  are  not  too  localized.  This  work  was 
reported  in  "Numerical  implementation  of  a  variational  method  for  electrical  impedance  tomogra¬ 
phy,"  by  R.  Kohn  and  A.  McKenney,  Inverse  Problems  b,  1990, 389-414. 

J.  Sylvester  explored  the  imaging  of  locally  anisotrophic  conductivities.  Here  there  is  a 
problem  of  non-uniqueness:  different  bodies  can  give  the  same  boundary  measurements  if  they 
are  related  -  "by  change  of  variables."  Sylvester  proved  the  converse,  in  two  dimensions:  two 
bodies  that  are  indistinguishable  by  boundary  measurements  must  be  related  in  this  way.  His 
work  is  reported  in  "An  anisotropic  inverse  boundary  value  problem,"  by  J.  Sylvester,  Comm. 
Pure  Appl.  Math.  43, 1990, 201-232. 

Other  Inverse  Problems 

C.S.  Morawetz  studied  two  numerical  problems  in  inverse  scattering  in  two  space  and  one 
time  dimensions  for  the  wave  equation.  The  basic  solutions  are:  in  free  space  a  plane  wave  for 
problem  I  and  a  fundamental  solution  for  problem  II.  Problem  I  is  the  classical  problem  of  deter¬ 
mining  a  potential  from  data  at  infinity  scattered  back  in  the  direction  from  which  the  wave  came. 
All  directions  of  the  initial  wave  then  give  enough  data.  Problem  n  is  to  determine  a  smooth 
sound  speed  form  data  collected  on  a  fixed  line  for  all  time.  This  problem  is  ill-posed  and  one 
can  only  ask  what  if  anything  can  be  determined  about  the  sound  speed. 

The  corresponding  direct  problems  are  analyzed  through  computed  examples.  In  Problem  I 
some  very  interesting  properties,  discovered  by  Ralston  and  Eskin,  of  the  scattered  field  interfere 
with  the  computation  of  the  inverse  and  it  is  impossible  to  make  a  very  logical  looking  algorithm 
converge.  In  Problem  n  for  the  direct  problem  a  back  scattered  field  was  obtained  showing  a 
wave  running  across  the  line.  The  first  part  of  this  work  was  reported  in  "Scattering  by  a  Potential 
Using  Hyperbolic  Methods",  by  A.  Bayliss,  Y.  Li  and  C.S.  Morawetz,  Mathematics  of  Computa¬ 
tion,  Vol.  52,  Number  186,  pp.  321-388,  (1989).  The  second  part  in  "On  2D  Inverse  Scattering", 
to  appear  in  The  Proceedings  of  the  International  Conference  on  the  Mathematical  and  Numeri¬ 
cal  Aspects  cfWcve  Propagation,  to  be  published  in  March  1991  by  SIAM. 
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Crystal  Growth 

Crystal  growth  computations  were  made  by  John  Strain  with  James  Sethian  (UC  Berkeley). 
They  wrote  and  debugged  a  code  for  solving  a  Stefan  problem  with  curvature-dependent  boun¬ 
dary  conditions.  The  code  uses  a  level  set  formulation  of  the  moving  boundary  problem,  with 
heat  potential  theory  to  reduce  the  heat  equation  to  the  moving  boundary.  A  fast  algorithm  is 
used  to  evaluate  the  heat  potential,  yielding  an  optimal  algorithm  which  is  capable  of  treating 
completely  arbitrary  geometries.  This  research  is  described  in  the  forthcoming  paper  "Crystal 
growth  and  Unstable  Solidification”. 

Strain  also  developed  a  new  family  of  fast  algorithms  for  classical  potential  theory,  based 
on  Ewald  summation.  These  algorithms  evaluate  classical  layer  or  volume  potentials,  or  discrete 
sums,  formed  with  the  Green  function  for  a  cube.  This  research  is  described  in  "Fast  Potential 
Theory  I:  Poisson  Solvers  on  a  Cube,"  Center  for  Pure  and  Applied  mathematics  (UC  Berkeley) 
Report  PAM-480,  December  1989,  and  submitted  to  Math.  Comp,  and  "Fast  Potential  H:  Layers 
Potentials  and  Discrete  Sums,"  Center  for  Pure  and  Applied  Mathematics  (UC  Berkeley)  Report 
PAM,  April  1990,  submitted  to  Jour.  Comput.  Phys. 

Further  work  has  been  on  quasistationary  crystal  growth  problems  using  integral  equations, 
level  sets  and  the  fast  algorithms  described  above  as  well  as  the  construction  of  a  fast  Laplace 
transform  based  on  Laguene  functions;  it  is  described  in  "A  Fast  Laplace  Transform  Based  on 
Laguerre  Functions,"  submitted  to  Math.  Comp. 

With  Leslie  Greengard,  he  constructed  a  fast  algorithm  for  evaluating  Gaussian  convolution 
sums,  "The  Fast  Gauss  Transform,"  to  appear  in  SIAM  J.  Sci.  Stat.  Comput.  This  work  was 
extended  to  evaluate  Gaussian  convolution  sums  with  variances  which  vary  from  point  to  point, 
"A  Fast  Gauss  Transform  with  Variable  Variances,"  submitted  to  SIAM  J.  Sci.  Stat.  Comput. 

Fast  Monte-Carlo  Methods: 

J.  Goodman,  with  A.  Sokal,  completed  a  study  with  theory  and  numerical  experiments  on 
their  multiscale  algorithm,  Multigrid  Monte-Carlo,  for  rapid  simulation  of  random  fields.They 
also  generalized  the  method  to  other  models,  such  as  the  XY  model,  with  significant  success. 
These  are  reported  in,  "Multi-Grid  Monte  Carlo  I.  Conceptual  Foundations”,  J.  Goodman  and  A. 
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D.  Sokal,  Phys.  Rev  A.  (1989)  and  "Multi-Grid  Monte  Carlo  II.  Two-Dimensional  XY  Model”, 
R.  Edwards  ,  J.  Goodman  and  A.  Sokal,  to  appear  in  Nuclear  Physics  B. 

Monte  Carlo  Methods  in  Phase  Transition 

Paula  A.  Whitlock  has  focussed  on  the  simulation  of  quantum  systems  using  Monte  Carlo 
methods.  The  distinct  areas  have  been  1)  the  simulation  of  helium  films  physisorbed  on  a  sub¬ 
strate;  2)  first  principles  calculation  of  two  and  three  body  interactions  and  3)  development  of 
Green’s  function  Monte  Carlo  method  employing  the  recently  developed  "shadow"  trial  wave 
function  for  quantum  many-body  systems.  All  these  areas  share  a  common  use  of  development 
of  quantum  Monte  Carlo  methods. 

The  ground-state  properties  of  liquid  and  solid  helium  are,  to  a  laige  extent,  well  under¬ 
stood.  However,  helium  systems  with  interfaces,  while  extensively  studied  experimentally,  are 
just  beginning  to  be  understood  by  theoreticians.  Whitlock  is  systematically  studying  helium 
adsorbed  on  graphite  using  Monte  Carlo  methods.  Helium  films  are  particularly  simple  and  their 
basic  properties  such  as  film  growth,  the  characteristics  of  the  films  two-dimensional  phases,  and 
their  melting-freezing  transitions  can  be  investigated  theoretically  and  compared  wtih  these 
phenomena  theoretically  at  the  microscopic  level,  as  well  as  calculating  from  microscopic  theory 
the  parameters  that  go  into  empirically-derived  theories. 

The  approach  to  the  study  of  helium  films  is  to  perform  variational  Monte  Carlo  and 
Green’s  function  Monte  Carlo  (GFMC)  simulations.  The  GFMC  method  was  first  introduced  by 
Kalos  and  was  further  refined  by  Kalos,  Levesque  and  Veriet  The  form  of  GFMC  used 
corresponds  to  resolvent  methods  in  many-body  perturbation  theory.  The  relevant  equation  is 
iterated  many  times,  it  converges  exponentially  and  exactly  to  the  lowest  state  of  the  Hamiltonian 
which  is  not  orthogonal  to  the  starting  wave  function.  The  iterations  are  carried  out  by  random 
walks  starting  with  an  optimized  trail  wave  function.  No  other  method  is  available  for  many- 
body  systems  which  can  evaluate  the  ground  state  and  orthogonal  excited  states  of  the  Schrod- 
inger  equation  exactly  within  the  statistical  sampling  errors.  The  results  of  the  calculation  are 
sets  of  particle  configurations  whose  probability  density  represent  the  ground  state  of  the  system. 
These  particle  configuration  are  used  to  form  expectation  values  to  compute  observables  of  the 
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system. 

A  variety  of  phenomena  associated  with  helium  films  have  been  studied.  The  equation  of 
state  can  be  computed  and  an  investigation  of  the  solid  phases  of  the  helium  monolayer  commen¬ 
surate  and  incommensurate  with  the  underlying  graphite  lattice  is  possible. 

Whitlock  also  worked  on  determining  the  properties  of  two  and  three  body  systems  through 
solving  the  Schrodinger  equation  exactly  using  Green’s  function  Monte  Carlo.  Whitlock’s 
approach  to  the  problem,  which  emphasizes  the  role  of  certain  interaction  energies  naturally 
divides  into  three  parts.  The  first  is  the  development  of  optimized  wave  functions  to  use  as  trial 
wave  functions  in  our  Monte  Carlo  simulations.  The  second  part  is  to  perform  variational  Monte 
Carlo  calculations  with  the  proposed  wavefunctions  and  the  last  part  is  to  perfom  the  Green’s 
function  Monte  Carlo  simulations.  So  far,  efforts  have  been  on  obtaining  the  two  and  three  body 
potentials  for  polarized  hydrogen  and  helium. 

Parallelization  of  algorithms 

This  work  was  conducted  mainly  under  the  guidance  of  J.  Demmel  and  in  the  first  year, 
Marsha  Berger.  It  involved  as  well  A.  Greenbaum,  A.  Mayo,  Z.  Bai,  A.  McKenney,  O.  Percus,  J. 
Schmidt. 

Much  of  the  work  revolved  around  the  production  of  the  LAPACK  linear  algebra  library  for 
parallel  computers.  This  is  due  to  be  released  shortly.  Its  goal  has  been  to  design  and  implement 
a  portable  linear  algebra  library  for  efficient  use  of  high-performance  supercomputers.  The  major 
methodology  for  making  the  algorithms  run  faster  was  to  restructure  them  to  perform  block 
matrix  operations  in  their  inner  loops.  These  functions  can  be  optimized  to  exploit  the  memory 
hierarchy  and  gain  large  speed  ups. 

Another  area  has  been  speeding  up  by  parallelization  well  known  codes  for  linear  systems 
arising  from  integral  equations  as  well  as  from  elliptic  partial  differential  functions.  The  first  pro¬ 
duces  dense  and  therefore  more  difficult  linear  systems  -  although  important,  too  difficult  to  have 
been  treated  much  in  the  past  even  though  they  are  well  conditioned.  Work  was  also  done 
analyzing  a  discrete-time  network  of  queues  under  various  circumstances.  This  is  an  important 
ingredient  in  making  a  parallel  system  work. 
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NumericaJ  Linear  Algebra 

James  Demmel  worked  on  algorithms  for  numerical  linear  algebra,  both  from  the  point  of 
view  of  improving  their  speed  on  supercomputers,  as  well  as  improving  their  accuracy.  This 
work  was  done  jointly  with  Zhaojun  Bai,  Anne  Greenbaum  and  Alan  McKenney.  They  have  been 
part  of  a  much  larger  team  including  researchers  from  U.  of  Tennessee  and  the  Numercial  Algo¬ 
rithms  Group,  Ltd.  This  is  the  work  described  above.  Other  results  were: 

"Jacobi’s  method  is  more  accurate  than  QR"  (by  J.  Demmel  and  K.  Veselic,  submitted  to 
SIAM  J.  Matrix  Analysis  and  Applications  ).  Shows  that  Jacobi’s  method  (with  a  modified  stop¬ 
ping  criterion)  computes  small  eigenvalues  of  symmetric  positive  definite  matrices  with  uni¬ 
formly  higher  relative  accuracy  than  QR,  divide  and  conquer,  tradition  bisection,  or  any  algo¬ 
rithm  which  first  involves  tridiagonalizing  the  matrix.  In  fact,  modulo  an  assumption  based  on 
extensive  numerical  tests,  we  show  that  Jacobi’s  method  is  optimally  accurate  in  the  following 
sense:  if  the  matrix  entries  have  small  relative  errors,  so  that  the  eigenvalues  have  certain  intrinsic 
uncertainties,  Jacobi  will  compute  them  with  nearly  this  accuracy.  In  other  words,  as  long  as  the 
initial  matrix  has  small  relative  errors  in  each  component,  even  using  infinite  precision  will  not 
improve  on  Jacobi  (modulo  factors  of  dimensionality).  Does  prove  bisection  is  optimal  in  this 
sense.  Also  shows,  the  eigenvectors  are  computed  more  accurately  by  Jacobi  (and  inverse  itera¬ 
tion)  than  previously  thought  possible.  Proves  similar  results  for  using  one-sided  Jacobi  for  the 
singular  value  decomposition  of  a  general  matrix.  Presents  a  version  of  Jacobi  with  the  property 
that  the  more  its  accuracy  exceeds  that  of  QR,  the  faster  it  converges. 

"The  Bidiagonal  Singular  Value  Decomposition  and  Hamiltonian  Mechanics"  (by  P.  Deift,  J. 
Demmel,  L.-C.  Li,  C.  Tomei,  accepted  by  SIAM  J.  Numerical  Analysis ).  Considers  the  problem 
of  computing  the  singular  value  decomposition  of  a  bidiagonal  matrix  B  .  This  problem  arises  in 
the  singular  value  decomposition  of  a  general  matrix,  and  in  the  eigenproblem  for  a  symmetric 
positive  definite  tridiagonal  matrix.  Shows  that  if  the  entries  of  B  are  known  with  high  relative 
accuracy,  the  singular  values  and  singular  vectors  of  B  will  be  determined  to  much  higher  accu¬ 
racy  than  the  standard  perturbation  theory  suggest  Also  shows  that  the  algorithm  in  [Demmel 
and  Kahan]  computes  the  singular  vectors  as  well  as  the  singular  values  to  this  accuracy.  Gives  a 
Hamiltonian  interpretation  of  the  algorithm  and  use  differential  equation  methods  to  prove  many 
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of  the  basic  facts.  The  Hamiltonian  approach  suggests  a  way  to  use  flows  to  predict  the  accumu¬ 
lation  of  error  in  other  eigenvalue  algorithms  as  well. 

"On  a  Block  Implementation  of  Hessenberg  Multishift  QR  Iteration"  (by  Z.  Bai,  J.  Demmel, 
appeared  in  International  J.  of  High-Speed  Computing,  v.  1,  n.  1, 1989).  The  usual  QR  algorithm 
for  finding  the  eigenvalues  of  a  Hessenberg  matrix  H  is  based  on  vector-vector  operations,  e.g. 
adding  a  multiple  of  one  row  to  another.  The  opportunities  for  parallelism  in  such  an  algorithm 
are  limited.  This  paper  describes  a  reorganizaton  of  the  QR  algorithm  to  peimit  either  matrix- 
vector  or  matrix-matrix  operations  to  be  performed,  both  of  which  yield  more  efficient  implemen¬ 
tations  on  vector  and  parallel  machine ,.  The  idea  is  to  chase  a  k  by  k  bulge  rather  than  a  1  by  1 
or  2  by  2  bulge  as  in  the  conventional  QR  algorithm.  Preliminary  numerical  experiments  on  the 
CONVEX  C-l  and  CYBER  205  vector  machines  are  reported. 

On  the  Conditioning  of  the  Nonsymmetric  Eigenproblem:  Theory  and  Software  by  Z.  Bai,  J. 
Demmel,  A.  McKenny.  This  report  reviews  the  theory  and  practical  estimation  of  condition 
numbers  of  the  nonsymmetric  eigenvalue  problem.  Also  provides  a  manual  for  using  LAPACK 
subroutines  STRSNA  and  STRSEN  to  estimate  condition  numbers  for  individual  eigenvlues  and 
eigenvectors,  multiple  (or  clustered)  eigenvalues,  and  invariant  subspaces. 

Iterative  methods 

In  addition  to  the  work  described  above  Arme  Greenbaum  has  dealt  with  several  different 
application  areas  and  has  focused  on  the  use  of  iterative  methods  to  solve  the  linear  systems  that 
rise  in  these  different  settings  and  on  the  use  of  parallel  computing  to  enable  more  rapid  solution. 
She  has  analyzed  various  preconditioners  for  the  linear  systems  arising  from  self-adjoint  elliptic 
partial  differential  equations  and  has  worked  on  practical  iterative  methods  for  the  solution  of 
nonsymmetric  linear  systems.  Much  of  the  work  has  dealt  with  iterative  solution  of  the  dense 
linear  systems  that  arise  from  integral  equations.  This  is  an  area  that  has  not  received  much 
attention  from  numerical  linear  algebraists  in  the  past  but  which  is  sure  to  prove  more  important 
in  the  future.  The  matrices  that  arise  here  are  dense  and  nonsymmetric  but  very  well-conditioned. 
This  means  that  iterative  methods  can  be  expected  to  converge  rapidly  and  the  major  part  of  the 
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work  will  be  in  applying  the  matrix  to  a  given  vector  at  each  step  of  the  iteration.  This  can  be 
done  rapidly  using  the  fast  multipole  method  of  Greengard  and  Rokhlin.  Greenbaum  wo  iced 
with  A.  Mayo  in  applying  these  solution  techniques  to  magnetics  problems  involving  Laplace’s 
or  the  biharmonic  equation  on  an  irregular  region.  Sh*  began  work  in  applying  such  techniques 
to  solve  a  free  boundary-value  problem  simulating  the  process  of  Ostwald  ripening.  She  has 
parallelized  the  adaptive  version  of  the  multi  pole  algorithm  on  the  NYU  Ultracomputer  proto¬ 
type.  Similar  techniques  are  being  applied  to  solve  scattering  problems  described  by  the 
Helmhotz  equation  on  an  exterior  region. 

Queuing  The  ~v 

In  problems  of  queuing  arising  from  parallelization  Ora  Pcrcus  worked  on 

i)  constructing  and  solving  carefully  chosen  models  in  order  to  illuminate  both  local  and 

global  structure  of  steady  state  operation  of  clock-regulated  queueing  networks  and  on 

ii)  designing  pseudorandom  number  generators  for  current  and  future  parallel  processors. 

Switching  networks  in  large  scale  digital  computers  are  characterized  by  synchronous 
pulsed  operation,  and  are  substantially  more  difficult  to  analyze  than  traditional  continuous  time 
queueing  systems.  The  control  parameter  in  such  a  network  is  the  probability  p  of  a  packet  enter¬ 
ing  a  channel  per  clock  cycle;  in  steady  state,  with  nonsaturating  queues,  p  is  constant  throughout 
the  network.  The  quantity  of  practical  interest  is  the  distribution  of  queue  length  throughout  the 
network,  since  a  queue  which  does  saturate  must  enter  another  mode  to  allow  continued  opera¬ 
tion.  It  is  important  to  have  a  concise  form  for  the  time  series  transformation  engendered  by  units 
of  queuing  networks;  the  combination  of  one  mixer  and  two  output  queues,  termed  a  2  x  2 
switch,  is  the  basic  unit  in  our  initial  investigation. 

Since  these  investigations  did  point  strongly  to  the  fact  that  the  near  independent  "light 
traffic"  assumption  was  a  surprisingly  decent  zeroth  approximation  up  to  half  of  the  saturation 
flow,  we  continued  by  investigating  the  distribution  of  the  queue  length  by  an  expansion  of  the 
time  series  in  a  channel  about  the  regime  of  independent  arrivals.  Results  for  the  second  stage  of 
such  a  network  that  we  obtained  recently  are  in  very  good  agreement  with  computer  simulations. 
This  points  *o  the  possiblity  that  the  i.i.d.  input  assumption  that  is  common  to  virtually  many 
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discrete  queueing  problems  can  be  extended  in  similar  fashion.  In  addition,  the  techniques  we 
developed  for  obtaining  the  queue  length  distribution  for  the  second  state  are  expected  to  give  us 
nontrivial  hints  as  to  how  to  work  at  higher  stages.  The  work  of  designing  pseudorandom  number 
generators  (pm’s)  for  current  and  future  parallel  processors  was  pointed  towards  generators  for 
MIMD  (multiple  instruction  stream/multiple  data  stream)  machines  such  as  the  Ultracomputer. 
Such  architectures  are  very  well  suited  for  Monte  Carlo  calculations,  since  potentially  indepen¬ 
dent  realization  or  "histories"  may  be  followed  on  each  processor.  They  are  much  more  suitable 
than  the  current  generation  of  "vector’  supercomputers  for  those  Monte  Carlo  calculations  with 
substantial  logical  complexity,  i.e.,  those  with  data-dependent  or  random  number  determined 
branches.  The  issues  addressed  are  applicable  to  many  parallel  machines,  including  existing 
vector-parallel  machines  (such  as  the  various  Crays,  multipipe  Cyber  205’s,  and  the  IBM  3094- 
400,  and  various  hypercube  machines).  The  emphasis  has  been  on  problems  associated  with  very 
large  numbers  of  processors. 

Issues  related  to  the  design  of  pseudorandom  number  generators  for  MIMD  parallel  proces¬ 
sors,  which  are  particularly  appropriate  for  Monte  Carlo  calculations  were  also  analyzed  to  ensure 
reproducibility  of  runs,  to  provide  very  long  sequences,  and  to  assure  an  adequate  degree  of 
independence  of  the  parallel  stream:,. 

Other  Topics 

The  streamline  diffusion  finite  element  method 

This  has  been  examined  by  Anders  Szepessy  and  is  a  higher  order  accurate  implicit  t>me 
stepping  scheme  for  hyperbolic  problems  on  general  meshes.  In  previous  work  he  had  proved 
convergence  of  these  approximations  for  a  general  scalar  conservation  law  in  several  space 
dimensions,  including  also  boundary  conditions.  It  is  assumed  that  the  nonlinear  discrete  equa¬ 
tions  on  each  time  step  are  solved  exactly. 

During  the  past  year,  Szepessy  analysed  and  modified  this  method  for  systems  of  conserva¬ 


tion  laws. 
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Since  the  equations  are  nonlinear  and  the  Unite  element  method  being  considered  is  impli¬ 
cit.  an  iterative  method  has  to  be  used  to  solve  the  discrete  (nonlinear)  equations  at  each  time 
step.  By  analyzing  the  viscous  profile,  resulting  from,  the  discretization,  Szepessy  proves  that  the 
above  convergence  result  for  scalar  problems  also  holds  when  the  nonlinear  discrete  equations, 
with  n  degrees  of  freedom,  are  solved  approximately  by  doing  a  few  iterations  (clog  n)  on  each 
time  step  with  a  simple  quasi-Newton  method.  To  obtain  this,  a  local  CFL-condition  has  to  be 
satisfied.  If  the  exact  solutions  of  the  conservation  law  is  smooth,  then  the  iterative  method  gives 
the  same  order  of  accuracy  as  solving  the  nonlinear  discrete  equations  exactly. 

Further,  by  diagonalizing  and  mollifying  the  artificial  viscosity  matrix  in  the  characteristic 
directions,  he  showed  that  viscous  profiles  also  exist  for  the  approximate  solutions  of  systems  of 
conservation  laws.  This  modification  improves  the  numerical  results  considerably. 

Vortex  Methods 

The  theory  of  vortex  methods  had  always  been  inapplicable  to  practical  computations  for 
two  reasons:  it  assumed  much  more  smoothing  than  is  actually  practical,  and  it  relied  on  Eulerian 
derivatives  of  the  Lagrangian  variable.  These  derivatives  grow  very  rapidly  even  for  tame  flows 
such  as  a  vortex  patch  and  lead  to  wildly  pessimistic  error  bounds.  Both  of  these  difficulties  have 
been  eliminated  in  joint  work  with  J.  Goodman.  First,  with  Cottet,  Hou,  and  Lowengrub,  conver¬ 
gence  of  the  vortex  method  with  no  smoothing  at  all  is  established  (the  point  vortex  method). 
These  convergence  theorems  agree  with  previous  numerical  tests  and  contradict  a  popular  body 
of  dogma.  Second,  with  Hou,  convergence  theorems  in  weak  nonns  that  depend  only  on  Eulerian 
quantities  are  established.  Numerical  experiments  indeed  show  that  errors  in  negative  norms 
grow  much  more  slowly  in  time  than  errors  in  strong  norms.  These  results  are  contained  in  "Con¬ 
vergence  of  the  Point  Vortex  Method",  J.  Goodman  and  T.  Hou,  to  appear  in  Comm.  Pure  Appl. 
Math,  "Convergence  of  a  Point  Vortex  Method  in  3D  with  Grid  Free  Stretching"  G.-H.  Cottet,  J. 
Goodman  and  T.  Hou,  to  appear  in  SIAM  J.  Num.  Anal,  and  J.  Goodman  and  T.  Hou  "New  Stabil¬ 
ity  Estimate  for  the  2-D  Vortex  Method"  submitted  to  Comm.  Pure  Appl.  Math. 
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Viscous  Shock  Waves 

Together  with  Zhouping  Xin,  Goodman  solved  the  "viscosity  method  problem"  for  piece- 
wise  smooth  solutions  of  systems  of  conservation  laws.  They  also  improved  the  stability  theory 
for  viscous  shock  waves  in  one  or  two  space  dimensions  through  a  number  of  technical  improve¬ 
ments:  the  use  of  Lp  norms  instead  of  L  2  and  the  derivation  of  the  "integrated  equations",  using  a 
"flux  variable  transformation  rather  than  by  integrating,  and  using  "interaction  weighted"  norms 
to  replace  the  complicated  space-time  estimates  of  the  earlier  theory.  See  "Viscous  limits  for 
Piecewise  Smooth  Solutions  to  Systems  of  Conservation  Laws",  in  preparation  and  "Remarks  on 
the  Stability  of  Viscous  Shock  Profiles",  preprint. 
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